# new from the version in gamlss.ts 
# the disSmo() functions disSmoA() and disSmoQ() can deal with NA in y
# I am not sure if this is the case for whitAR1() and whitAR2()
#------------------------------------------------------------
# AUTHORS 
# Mikis Stasinopoulos, Bob Rigby, Paul Eilers, Vlasis Voudouris and Majid  Djennad
# this file contains  functions to fit discrete smoother for time series data 
# It has following functions 

# vi)  WhitAR1 
# 1) disSmoA() : for "Alternate"" fitting when the smoothing parameters is estimated
# 2) disSmoQ() : for "Q function"" fitting when the smoothing parameters is estimated
#                 It produce also standard errors for the parameters
#                 functions disSmoA and disSmoQ  should give identical results
#                 Both functions are using  the new spam 0.24 package (Not in CRAN yet)
# 3) disSmo()   : which combines the above two functions (now disSmoQ is the default ) 
# 4) methods for an dissmo onject
#      i)   fitted()
#      ii)  coef()
#      iii) residuals()
#       iv) AIC  
#       v)  deviance
#       vi) plot
#     vii)  print()
#     viii) summary()
# 5) AR models
#           i) whitAR1() for AR(1) discrete Smoothing models
#          ii) whitAR2() for AR(2) discrete Smoothing models
#
# 7) eigenValues() : the eigen values for the P matrix of a RW1 
# 8 ) HessianPB() : calculate the hessian used after "nlminb" maximazion
#
#
#---------------------------------------------------------------
# IMPORTANT
# ** NOTE AT THE MOMENT CAN NOT RUN WITHOUT spam 0.24
#---------------------------------------------------------------
# The method in disSmoA "Alternating" between sigma and sigma_b 
# note that this method requires and estimates of the degrees of freedom
# The  degrees of freedom are the trace of the smoothing matrix
#       (W+lambdaG)^(-1)W 
# which is difficult to calculate for large N 
# therefore an approximation is used here.
# this approximation is very accurate
# for unweighted observations but can result to discrepancies if the data are weighted
# (The method compared well with the Q method (at least for unweihted data)
# where the df are esimated using formulae)
#--------------------------------------------------------------- 
# both methods create an "disSmo" object which has the following methods
# print(), fitted(), resid(), AIC(), deviance(), plot(), coef()
# Not that wp() and dtop() can also be used
#---------------------------------------------------------------
# NOTE : the limit of the functions disSmo() if lambda -> Infinity 
# is order d-1 polynomial
# for example
# order 1 -> constant
# order 2 -> linear
# order 3 -> quadratic 
# order 4 -> cubic
# Interporlation is done by an 2d-1 polynomial
# Extrapolation by a d-1 polynomial
#-------------------------------------------------------------------------------------
# questions
# do we allow order=0? (random effects)
#-------------------------------------------------------------------------------------
#  TO DO
# i)    predict() and  forecast() functions 
# ii)   missing values : missing values needs y and weights to be set to zero: this is fixed for
#       for the disSmo functions
# iii)  Standard Errors for the sigmas (and AR's'): this is fixed, they are calculated now for 
#         function which use the Q function 
# vi)   Q function does not always converge in the same value as the A with the same starting values
#        Note that in the disSmoQ function they exist now 3-alternatives ways for optimation
#        optim(), nlminb() and genoud(). The first and third provides hessian.  
#        The  nlminb (I think does not) so we use the function HessianPB() which cames from the nlme 
#        package
#        we need for expreriance with the genoud function so we know how to tune it
# v)    Correct for  weights in the weighted version
# vi)   version of disSmo() cycling in the end as in Rue and Held (2005) : Majid hav done it we need 
#       to check it 
# v)    standard errors for the fitted functions ??
# vi)   whitAR2 and whitAR3 need checking  
# vii)  We need RD files for all the functions 
#-------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------
# 1st  function
#-------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------
# METHOD ALTER : alternate ML estimation
disSmoA <- function(y, 
             weights = rep(1,length(y)),  
              lambda = NULL,
               order = 1,
               start = 10,
                plot = FALSE)                    
 {
 # local functions
 #------------------------------- 
      regpen <- function(y, weights, lambda)  
       {
               fv <- solve(W+lambda*G,  weights*y)
               if (order==0)
                 { # if order zero we can get them easy  
                    edf <- sum(weights/(weights+lambda))    
                 }
                else
                 { # here is more difficult to get dfs
                    if (N<=500) edf <- sum(diag(solve(W+lambda*G,  W)))
                    else  # otherwise get a sample of N1 
                     {    # and try to estimate it see   P. Eilers (2003) method (A Perfect Smoother)
                                     II <- seq(1,N, ceiling( N/(log(N)*50))) # get a sample 
                                     N1 <- length(II)                        # of length N1
                                     z2 <- rep(0,N)                          # create two dummy variable
                                     z1 <- rep(0,N1)
                    z2[as.integer(N/2)] <- 1                                 # put 1 in the middle
                   z1[as.integer(N1/2)] <- 1                                  
                                lambda1 <- lambda #((N1/N)^(2*order))*(lambda) # 
                                     #yy <- y[II]                             # get sample data
                                     ww <- weights[II] 
                                     WW <- diag.spam(x=ww)
                                     EE <- diag.spam(N1)     
                                     DD <- diff(EE, diff = order)   
                                     GG <- t(DD)%*%DD
                                   edf1 <- sum(H1<-diag(solve(WW+lambda1*GG, WW))) # get edf from the sample
                                    edf <- edf1 +(N-N1)* mean(H1[c(as.integer((N1/2)-100):as.integer((N1/2)+100))])
                  
                     } 
                  }      
             fit <- list(fv=fv, edf=edf)
            fit 
        }
 #-------------------------------	 
      require(spam)
      require(gamlss)
        scall <- deparse(sys.call())
       natrue <- FALSE
           if (any(is.na(y)))
           {            
                        natrue <- TRUE
                            yy <- y # we need this for plotting
             weights[is.na(y)] <- 0
                   y[is.na(y)] <- 0 
           }
           if (any(weights==0))
           {
                       natrue <- TRUE
                           yy <- ifelse(weights==0, 0, y) 
           }
            N <-  length(y)
            W <- diag.spam(x=weights)
            E <- diag.spam(N)     
            D <- diff(E, diff = order)   
            #D <- if (order==0) E
            #     else  -1*E[-N, ] + E[-1, ]
                 #as.spam(diff(E, differences=order))
            G <- t(D)%*%D 
            sig2 <- tau2 <- 0
            
         if (is.null(lambda)) 
           {
                    lambda <- start
                for (it in 1:200) 
                {
                          fit <- regpen(y, weights, lambda)
                       gamma. <- diff(fit$fv, differences=order)
                           fv <- fit$fv          
                         sig2 <- sum(weights * (y - fv) ^ 2) / (N - fit$edf)
                         tau2 <- sum(gamma.^2)/(fit$edf - order)
                   lambda.old <- lambda
                       lambda <- sig2/tau2
                  #cat("lambda", lambda, "sig2", sig2, "tau2", tau2,"\n")
                  if (lambda < 1e-07) lambda <- 1e-07
                  if (lambda > 1e+07) lambda <- 1e+07
                  if (abs(lambda - lambda.old) < 1e-07 || lambda > 1e+10) break
                }
              }            
           fit <- regpen(y, weights, lambda) 
          sig2 <- sum(weights * (y - fit$fv) ^ 2) / (N - fit$edf)
          tau2 <- sum(diff(fit$fv, differences=order)^2)/(fit$edf - order)   
           par <- c(sig2, tau2)
    names(par) <- c("sige^2", "sigb^2")
# saving staff 
           fit <- list(fitted = fit$fv, 
                           df = fit$edf, 
                       lambda = lambda, 
                        order = order, 
                         sig2 = sig2,
                        sigma = sqrt(sig2),
                         tau2 = tau2,
                        par = list(par=par, se=c(NA, NA)),
                            y = y,
                      weights = weights,
                            N = N,
                         call = scall,
                          rss = sum(weights*(y-fit$fv)^2),
                          aic = sum(weights*(-2*dNO(y, mu=fit$fv, sigma=sqrt(sig2), log=TRUE)))+2*(fit$edf+1) , 
                          sbc = sum(weights*(-2*dNO(y, mu=fit$fv, sigma=sqrt(sig2), log=TRUE)))+log(N)*(fit$edf+1),
                    deviance = sum(weights*(-2*dNO(y, mu=fit$fv, sigma=sqrt(sig2), log=TRUE))))
          class(fit) <- "disSmo"
           if (plot) 
           {
            if (natrue) plot(yy, col = gray(0.7), type = 'l')
            else plot(y, col = gray(0.7), type = 'l') 
           lines(fitted(fit)~as.numeric(time(y)), col = 'red', lwd = 2)
            }
           return(fit)     
}
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
# 2 function 
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------- 
# disSmoQ
# method Q function
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------      
# disctere smoothing using the Q function
# that is the evaluation of the smoothing parameters is done
# by maximising the Q-function of Pawitan 2002
# which is a profile likelihood for the sigma and tau parameters 
# it uses the spam package version 0.24 which is NOT in the CRAN yet
# this was based on the function disSmo3() of Mikis Stasinopoulos
#---------------------------------------------------------------------------------------
disSmoQ <- function(y, 
                weights = rep(1,length(y)), # for weighted observations 
                  order = 1,
                  start = 10,
                   plot = FALSE,
             optim.proc = c("optim", "nlminb", "genoud"))                       
 {          
  require(spam); require(gamlss)
        scall <- deparse(sys.call())
      optim.p <- match.arg(optim.proc)
        natrue <- FALSE 
            if (any(is.na(y)))
           {
                       natrue <- TRUE
                            yy <- y # we need this for plotting
             weights[is.na(y)] <- 0
                   y[is.na(y)] <- 0
           }
       if (any(weights==0))
           {
             natrue <- TRUE
             yy <- ifelse(weights==0, 0, y) 
           }
      lambda  <- start 
            N <- length(y)
            W <- diag.spam(x=weights)
            E <- diag.spam(N)
            D <- diff(E, diff = order)   
            G <- t(D)%*%D 
            #       BD <- bdiag(W+lambda*G) 
           fv <- solve(W+lambda*G,  weights*y)        
       params <- c(sige <- sum(weights*(y-fv)^2)/N, sigb<-sum((D%*%fv)^2)/N)        
    Qf<- function(par) 
     { # maybe we reparametrice to work on the log scale of the parameters 
      # cat("par", par, "\n")
       lambda <- par[1]/par[2]
           fv <- solve(W+lambda*G,  weights*y)   
            b <- diff(fv, differences=order)
           D3 <- determinant((1/par[1])*W+(1/par[2])*G)$modulus
            f <- -(N/2)*log(2*pi*par[1]) -sum(weights*(y-fv)^2)/(2*par[1])-((N-order)/2)*log(2*pi*par[2])-sum(b^2)/(2*par[2]) -.5*D3+(N/2)*log(2*pi)
           -f
     }

 switch(optim.p, "nlminb"={# this do not allow se's'
            out <- nlminb(start = params, objective = Qf, lower = c(1e-10, 1e-10),  upper = c(Inf, Inf))
    out$hessian <- HessianPB(pars=out$par, fun=Qf)$Hessian
     value.of.Q <- out$objective
                         },
                 "optim"={
            out<- optim(par= params, fn = Qf, lower = c(1e-10, 1e-10), upper = c(Inf, Inf),method= "L-BFGS-B",
                   hessian=TRUE)
                  # geting the hessian using nlme    fdHess(out$par, Qf)
            value.of.Q <- out$value
                           },
                  "genoud"={
                    require(rgenoud)
                    B <- matrix(c(1e-10,1e-10, 1e10, 1e10), ncol=2)
                    # spead is INCREASED  by reducing  pop.size and increasing BFGSburnin
              out <- genoud( Qf, nvar=2, max=FALSE, pop.size=20,  Domains=B, boundary.enforcement=2, gradient.check=TRUE, hessian=TRUE,BFGSburnin=20, max.generations=50,print.level=0)
            # out <- genoud( Qf, nvar=2, max=FALSE, pop.size=10,  Domains=B, boundary.enforcement=2, hessian=TRUE, BFGSburnin=50, P1=0, P2=50, P3=50, P4=0, P5=0, P6=0, P7=0, P8=0, P9=0)
            #out <- genoud( Qf, nvar=2, max=FALSE, pop.size=5000, Domains=B, boundary.enforcement=2, hessian=TRUE)
                    value.of.Q <- out$value   
                    })    
  #   out1<- optim(par= params,            fn = Qf, lower = c(1e-10, 1e-10),  upper = c(Inf, Inf),   method= "L-BFGS-B")
 #  library(rgenoud)
 #  B <- matrix(c(1e-10,1e-10, Inf, Inf), ncol=2)  
 # B <- matrix(c(0.001,0.001, 10000, 10000), ncol=2)
 #    genoud( Qf, nvar=2, max=FALSE, pop.size=3000, Domains=B)
      sige <- out$par[1]
      sigb <- out$par[2]
        se <- if (is.null(out$hessian)) rep(NULL,2) 
              else sqrt(diag(solve(out$hessian)))
names(out$par) <- c("sige^2", "sigb^2")
    lambda <- sige/sigb
        fv <- solve(W+lambda*G,  weights*y)   
         b <- diff(fv, differences=order)
       tr1 <- order + sum(b^2)/(sigb)
       tr2 <- N-(sum(weights*(y-fv)^2))/(sige)
    if (plot) {if (natrue) plot(yy, col = gray(0.7), type = 'l')
            else plot(y, col = gray(0.7), type = 'l')
               lines(fv~as.numeric(time(y)), col = 'red', lwd = 2)}
    # get the output 
           fit <- list(fitted = fv, 
                           df = tr1,
                           df2= tr2,
                       lambda = lambda, 
                        order = order, 
                         sig2 =  sige,
                        sigma = sqrt(sige),
                         tau2 = sigb,
                   value.of.Q = value.of.Q,
                          par = list(par=out$par,se=se), 
                            y = y,
                      weights = weights,
                            N = N,
                         call = scall,
                          rss = sum(weights*(y-fv)^2),
                          aic = sum(weights*(-2*dNO(y, mu=fv, sigma=sqrt(sige), log=TRUE)))+2*(tr1) , 
                          sbc = sum(weights*(-2*dNO(y, mu=fv, sigma=sqrt(sige), log=TRUE)))+log(N)*(tr1),
                     deviance = sum(weights*(-2*dNO(y, mu=fv, sigma=sqrt(sige), log=TRUE))))
          class(fit) <- "disSmo"
  fit
  }
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
# 3 function
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
# disSmo 
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
disSmo<- function(y, 
              method = c( "Qfun", "Alter"),
                   ...)
{
  method <- match.arg(method)
   scall <- deparse(sys.call())
     fit <- switch(method, "Alter"= disSmoA(y,...),
                    "Qfun" = disSmoQ(y,...)) 
fit$call <- scall
 fit
}
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
# functions 4 disSmo methods
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
# methods for disSmo objects                    
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------                    
fitted.disSmo<-function(object,...) 
{
object$fitted
}
#------------------------------------------
#------------------------------------------
coef.disSmo<-function(object,...) 
{
coef <- object$par[["par"]]
attr(coef, "se") <-  object$par[["se"]]
coef
}
  
#coef.whitAR1<-function(object,...) 
#{
#object$par
#}
#------------------------------------------
#------------------------------------------
residuals.disSmo<-function(object,...) 
{
res<-qNO(object$weights*pNO(object$y, mu=object$fitted, sigma=sqrt(object$sig2)))
res <- ifelse(res==-Inf, NA, res)
res <- ifelse(res==Inf, NA, res)
res[!is.na(res)]
}
#-------------------------------------------
#------------------------------------------
AIC.disSmo <- function(object, ...,k=2)
{
 val <- if (is(object, "disSmo")) 
            object$deviance + (object$df+1) * k
        else stop(paste("this is not a disSmo object"))
val
}
#------------------------------------------
#------------------------------------------
deviance.disSmo<-function(object,...) 
{
object$deviance
}
#------------------------------------------
#------------------------------------------
plot.disSmo <- function (x, xvar=NULL, parameters=NULL, ts=FALSE, summaries=TRUE, ...) 
{
    residx <- resid(x) # get the residuals 
         w <- x$weights
       x$x <- if (is(x$y , "ts"))  as.numeric(time(x$y))
              else                 1:length(x$y)
    xlabel <- if(!missing(xvar)) deparse(substitute(xvar)) else deparse(substitute(index))
## plotting parameters
    if(is.null(parameters))
          op <- par(mfrow=c(2,2), mar=par("mar")+c(0,1,0,0), col.axis="blue4", col.main="blue4", col.lab="blue4",  col="darkgreen", bg="beige" )
    else  op <- parameters
## now the two top  figures 
## if time series plot acf and pacf  
    if(identical(ts, TRUE))
     {  # get the acf and pacf
     acf.new<-acf(residx,plot=FALSE)
     plot(acf.new,xlim=c(2,length(acf.new$acf)),ylim=range(acf.new$acf[-1]))   # ms Tuesday, August 19, 2003 at 11:04
     pacf(residx)
     }
     else 
     {# otherwise 
     ## I am assuming that is x$noObs!=x$N then we have weights (with frequencies)
     if (length(residx)==x$N)
        {
         fittedvalues <- if(is.null(fitted(x))) fitted(x,"sigma") else fitted(x) # MS Wednesday, September 10, 2003 at 21:20
         ## whether index or x-variable
         if(is.null(xvar))     xvar <- seq(1,length(residx),1) # MS
        }
        else
        { # if weights
         fittedvalues <- rep( if(is.null(fitted(x))) fitted(x,"sigma") else fitted(x), w)
          xvar <- if(is.null(xvar))  seq(1,length(residx),1) else rep(xvar,w)
        } 
    # top left
    plot(fittedvalues , residx,
         xlab = "Fitted Values",  
         ylab = "Quantile Residuals", 
         main = "Against Fitted Values",
         frame.plot = TRUE) 
    # top right  
    plot(xvar, residx, 
         ylab = "Quantile Residuals",
         xlab = xlabel, 
         main = paste("Against ", xlabel), 
         frame.plot = TRUE) #  points(par(col="blue4"))
     }    
    plot(density(residx), 
         xlab = "Quantile. Residuals", 
         ylab = "Density", 
         main = "Density Estimate",
         frame.plot = TRUE, 
         col="black", 
         lwd=0.4 ) #col="deepskyblue4", col="darkgreen", 
         rug(residx, col="red")
 
    qqnorm(residx, main = "Normal Q-Q Plot",
            xlab = "Theoretical Quantiles",
            ylab = "Sample Quantiles", 
            plot.it = TRUE, 
            frame.plot = TRUE, 
            col="darkgreen")
       lines(as.numeric(residx), as.numeric(residx), col="red" , lwd=.4, cex=.4 )
 
     if ( identical(summaries, TRUE))
               { 
                     qq <- as.data.frame(qqnorm(residx, plot = FALSE))
               Filliben <- cor(qq$y,qq$x)
                    # mr <- as.matrix(residx)
                    m.1 <- mean(residx)
                    m.2 <- var(residx) # cov.wt(mr,w)$cov
                  n.obs <- sum(w) 
                    m.3 <- sum((residx-m.1)**3)/n.obs 
                    m.4 <- sum((residx-m.1)**4)/n.obs 
                    b.1 <- m.3^2/m.2^3
                sqrtb.1 <- sign(m.3)*sqrt(abs(b.1))
                    b.2 <- m.4/m.2^2 
                     cat("*******************************************************************")
                     cat("\n")
                     if (identical(x$type,"Continuous")) 
                         {cat("\t","     Summary of the Quantile Residuals")}
                     else{cat("\t","Summary of the Randomised Quantile Residuals")}    
                     cat("\n")
                     cat("                           mean   = ", m.1, "\n")
                     cat("                       variance   = ", m.2, "\n")
                     cat("               coef. of skewness  = ", sqrtb.1, "\n")
                     cat("               coef. of kurtosis  = ", b.2, "\n")
                     cat("Filliben correlation coefficient  = ", Filliben, "\n")
                     cat("*******************************************************************")
                     cat("\n")

               }
    par(op)
}
#------------------------------------------
#------------------------------------------
print.disSmo  <- function (x, digits = max(3, getOption("digits") - 3), ...) 
{   
    cat("\nDiscrete Smoothing fit")
    cat("Fitting method:", deparse(x$method), "\n")
    cat("\nCall: ", deparse(x$call), "\n", fill = TRUE)
    cat("Coefficients")
    co <- coef(x)
    cat("\n  ", names(co), "\n")
    cc <-simplify2array(co, higher=TRUE)
    cat(cc, " \n")
    cat("\n Degrees of Freedom for the fit:", x$df, "Residual Deg. of Freedom  ", 
        x$N-x$df, "\n")
    cat("Global Deviance:    ", format(signif(x$deviance)), 
        "\n            AIC:    ", format(signif(x$aic)), "\n            SBC:    ", 
        format(signif(x$sbc)), "\n")
    invisible(x)
}
   
summary.disSmo  <- function (object, digits = max(3, getOption("digits") - 3), ...) 
{   
    cat("\nDiscrete Smoothing fit")
    cat("Fitting method:", deparse(object$method), "\n")
    cat("\nCall: ", deparse(object$call), "\n", fill = TRUE)
            coef  <- coef(object)
          se.coef <- attributes(coef)$se
             tval <- coef/se.coef
          matcoef <- cbind(coef, se.coef, tval, 2*(1-pnorm(abs(tval))))
dimnames(matcoef) <- list(names(tval), c(" Estimate", " Std. Error", " t value", "Pr(>|t|)"))
cat("\nCoefficient(s):\n")
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)

   
    cat("\n Degrees of Freedom for the fit:", object$df, "Residual Deg. of Freedom  ", 
        object$N-object$df, "\n")
    cat("Global Deviance:    ", format(signif(object$deviance)), 
        "\n            AIC:    ", format(signif(object$aic)), "\n            SBC:    ", 
        format(signif(object$sbc)), "\n")
    invisible(object)
}
# methods are finish here    
#---------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
# 5 function AR  
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#18.07.11 AR(1) Majid function 
whitAR1 <- function(y, 
                   weights = rep(1,length(y)), # for weighted observations 
                     order = 1,
                     start = 10,
                      init = 0.5,
                      optim.proc=c("optim", "nlminb", "genoud"),
                      plot = FALSE)                       
 {          
 
  require(spam); require(gamlss)
          scall <- deparse(sys.call())
        optim.p <- match.arg(optim.proc)
         natrue <- FALSE 
            if (any(is.na(y)))
           {
                       natrue <- TRUE
                            yy <- y # we need this for plotting
             weights[is.na(y)] <- 0
                   y[is.na(y)] <- 0
           }
       if (any(weights==0))
           {
             natrue <- TRUE
             yy <- ifelse(weights==0, 0, y) 
           }   
  lambda  <- start 
           ro <- init
            N <- length(y)
            W <- diag.spam(x=weights)
            E <- diag.spam(N)
            D <- if (order==0) E
                 else  -ro*E[-N, ] + E[-1, ]          
            G <- t(D)%*%D  
           fv <- solve(W+lambda*G,   weights*y)        
          params <- c(sige <- sum( weights*(y-fv)^2)/N, sigb<-sum((D%*%fv)^2)/N, ro=ro) 
          names(params) <- c("sige^2", "sigb^2", "phi")
  #------------------------  
  Qf<- function(par)
      {  
        lambda <- par[1]/par[2]
             D <-  -par[3]*E[-N, ] + E[-1, ]
        #cat(par, "\n")
             G <-   t(D)%*%D
             fv <- solve(W+lambda*G,   weights*y)
             b <-   D %*% fv
             D3 <-  determinant((1/par[1])*W+(1/par[2])*G)$modulus
            f <- -(N/2)*log(2*pi*par[1]) - sum( weights*(y-fv)^2)/(2*par[1]) - ((N-order)/2)*log(2*pi*par[2]) - sum(b^2)/(2*par[2]) -.5*D3 + (N/2)*log(2*pi)
           -f
      }
 #------------------------
      switch(optim.p, "nlminb"={# this do not allow se's'
            out <- nlminb(start = params, objective = Qf, lower = c(1e-10, 1e-10, -0.999999999),  upper = c(Inf, Inf, .999999999))
            value.of.Q <- out$objective
            out$hessian <- HessianPB(pars=out$par, fun=Qf)$Hessian
                         },
                 "optim"={
            out<- optim(par= params, fn = Qf, lower = c(1e-10, 1e-10, -.999999999), upper = c(Inf, Inf, .999999999), method= "L-BFGS-B",
                   hessian=TRUE)
            value.of.Q <- out$value
                           },
                  "genoud"={
                    require(rgenoud)
                    B <- matrix(c(1e-10, 1e-10, -0.9999999, 1e10, 1e10, 0.9999999), ncol=2)
             out <- genoud( Qf, nvar=3, max=FALSE, pop.size=20,  Domains=B, boundary.enforcement=2, 
                            gradient.check=TRUE, hessian=TRUE,BFGSburnin=20, max.generations=50,
                            print.level=0)
             value.of.Q <- out$value   
                    })
      sige <- out$par[1]
      sigb <- out$par[2]
        phi <- out$par[3] 
        se <- if (is.null(out$hessian)) rep(NULL,2) 
              else sqrt(diag(solve(out$hessian)))
    lambda <- sige/sigb
         D <- -phi*E[-N, ] + E[-1, ]
         G <-   t(D)%*%D
        fv <- solve(W+lambda*G,   weights*y)
         b <-   D %*% fv 
       tr1 <- order + sum(b^2)/(sigb)
       tr2 <- N-(sum(weights*(y-fv)^2))/(sige)
    if (plot) {if (natrue) plot(yy, col = gray(0.7), type = 'l')
            else plot(y, col = gray(0.7), type = 'l')
               lines(fv~as.numeric(time(y)), col = 'red', lwd = 2)}
    # get the output 
           fit <- list(fitted = fv, 
                           df = tr1,
                           df2= tr2,
                       lambda = lambda,
                          par = list(par=out$par, se=se),
                           ro = ro,
                        order = order, 
                         sig2 = sige,
                        sigma = sqrt(sige),
                         tau2 = sigb,
                            y = y,
                      weights =  weights,
                            N = N,
                         call = scall,
                          rss = sum(weights*(y-fv)^2),
                          aic = sum(weights*(-2*dNO(y, mu=fv, sigma=sqrt(sige), log=TRUE)))+2*(tr1) , 
                          sbc = sum(weights*(-2*dNO(y, mu=fv, sigma=sqrt(sige), log=TRUE)))+log(N)*(tr1),
                     deviance = sum(weights*(-2*dNO(y, mu=fv, sigma=sqrt(sige), log=TRUE))))
       class(fit) <- c( "whitAR1", "disSmo")        
    fit
}
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#12.08.11 AR(2) Majid function 
# extending whittaker to AR(2)
whitAR2 <- function(y, 
                   weights = rep(1,length(y)), # for weighted observations 
                     order = 2,
                     start = 10,
                   ro.init = 0.5,
                  phi.init = 0.1,
                      plot = FALSE)                       
 {          
  require(spam); require(gamlss)
        scall <- deparse(sys.call())
#---------------------------------------------------------------------------
   natrue <- FALSE 
            if (any(is.na(y)))
           {
                       natrue <- TRUE
                            yy <- y # we need this for plotting
             weights[is.na(y)] <- 0
                   y[is.na(y)] <- 0
           }
       if (any(weights==0))
           {
             natrue <- TRUE
             yy <- ifelse(weights==0, 0, y) 
           }   
#---------------------------------------------------------------------------  
      lambda  <- start 
           ro <- ro.init
          phi <- phi.init
            N <- length(y)
            W <- diag.spam(x=weights)
            E <- diag.spam(N)
            D <- if (order==0) E
                 else E[-c(1,2), ] - ro*E[-c(1,N),  ] - phi*E[-c(N-1,N), ]  
            G <- t(D)%*%D  
           fv <- solve(W+lambda*G,   weights*y)        
         par <- c(sige <- sum( weights*(y-fv)^2)/N, sigb<-sum((D%*%fv)^2)/N, ro=ro, phi=phi) 
  names(par) <- c("sige^2", "sigb^2", "phi1", "phi2")
    Qf<- function(par)
      {
        lambda <- par[1]/par[2]
             D <-  E[-c(1,2), ] - par[3]*E[-c(1,N),  ] - par[4]*E[-c(N-1,N), ]
             G <-   t(D)%*%D
             fv <- solve(W+lambda*G,   weights*y)
             b <-   D %*% fv
             D3 <-  determinant((1/par[1])*W+(1/par[2])*G)$modulus
            f <- -(N/2)*log(2*pi*par[1]) - sum( weights*(y-fv)^2)/(2*par[1]) - ((N-order)/2)*log(2*pi*par[2]) - sum(b^2)/(2*par[2]) -.5*D3 + (N/2)*log(2*pi)
           -f
      }

    out<-nlminb(start = par, objective = Qf, lower = c(0.00001,0.000001, -100, -100),  upper = c(Inf, Inf, 100, 100))
      
      sige <- out$par[1]
      sigb <- out$par[2]
        ro <- out$par[3]
       phi <- out$par[4]
    lambda <- sige/sigb
         D <-  E[-c(1,2), ] - ro*E[-c(1,N),  ] - phi*E[-c(N-1,N), ]
         G <-   t(D)%*%D
        fv <- solve(W+lambda*G,   weights*y)
         b <-   D %*% fv 
       tr1 <- order + sum(b^2)/(sigb)
       tr2 <- N-(sum(weights*(y-fv)^2))/(sige)
    if (plot) {if (natrue) plot(yy, col = gray(0.7), type = 'l')
            else plot(y, col = gray(0.7), type = 'l')
               lines(fv~as.numeric(time(y)), col = 'red', lwd = 2)}
    # get the output 
           fit <- list(fitted = fv, 
                           df = tr1,
                           df2= tr2,
                       lambda = lambda,
                          par = list(par=out$par),
                           ro = ro,
                          phi = phi,
                        order = order, 
                         sig2 =  sige,
                        sigma = sqrt(sige),
                         tau2 = sigb,
                            y = y,
                      weights =  weights,
                            N = N,
                         call = scall,
                          rss = sum(weights*(y-fv)^2),
                          aic = sum(weights*(-2*dNO(y, mu=fv, sigma=sqrt(sige), log=TRUE)))+2*(tr1) , 
                          sbc = sum(weights*(-2*dNO(y, mu=fv, sigma=sqrt(sige), log=TRUE)))+log(N)*(tr1),
                     deviance = sum(weights*(-2*dNO(y, mu=fv, sigma=sqrt(sige), log=TRUE))))
       class(fit) <- c( "whitAR", "whitAR1", "disSmo")
          
    fit
  
   
}

# how to find the eigen values fot the RW1 G matrix
eigenValues<-function(N)# form Rue and Held (2005) page:96
{
  i<-1:N
  Lambda <- 2-2*cos(pi*((i-1)/N))
  Lambda
}
#-------------------------------------------------------------------
#-------------------------------------------------------------------
#-------------------------------------------------------------------
# this function from nlme of Pinheiro and Bates 
 HessianPB<-function (pars, fun, ..., .relStep = (.Machine$double.eps)^(1/3), 
    minAbsPar = 0) 
{
    pars <- as.numeric(pars)
    npar <- length(pars)
    incr <- ifelse(abs(pars) <= minAbsPar, minAbsPar * .relStep, 
        abs(pars) * .relStep)
    baseInd <- diag(npar)
    frac <- c(1, incr, incr^2)
    cols <- list(0, baseInd, -baseInd)
    for (i in seq_along(pars)[-npar]) {
        cols <- c(cols, list(baseInd[, i] + baseInd[, -(1:i)]))
        frac <- c(frac, incr[i] * incr[-(1:i)])
    }
    indMat <- do.call("cbind", cols)
    shifted <- pars + incr * indMat
    indMat <- t(indMat)
    Xcols <- list(1, indMat, indMat^2)
    for (i in seq_along(pars)[-npar]) {
        Xcols <- c(Xcols, list(indMat[, i] * indMat[, -(1:i)]))
    }
    coefs <- solve(do.call("cbind", Xcols), apply(shifted, 2, 
        fun, ...))/frac
    Hess <- diag(coefs[1 + npar + seq_along(pars)], ncol = npar)
    Hess[row(Hess) > col(Hess)] <- coefs[-(1:(1 + 2 * npar))]
    list(mean = coefs[1], gradient = coefs[1 + seq_along(pars)], 
    
        Hessian = (Hess + t(Hess)))
}
#------------------------------------------------------------------------